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Abstract 


The objective of this investigation has been to develop an algo- 
rithm (or algorithms) for the improvement of the accuracy and ef- 
ficiency of the computer fluid dynamics (CFD) models to study the 
fundamental physics of combustion chamber flows, which are neces- 
sary ultimately for the design of propulsion systems such as SSME 
and STME. During this three year study (May 19, 1989 - Mary 18, 

1992), a unique algorithm was developed for all speed flows. This 
newly developed algorithm basically consists of two pressure-based 
algorithms (i.e. PISOC and MFICE). This PISOC is a non-iterative 
scheme and the FICE is an iterative scheme where PISOC has the 
characteristic advantages on low and high speed flows and the modi- 
fied FICE has shown its efficiency and accuracy to compute the flows 
in the tr ans onic region. A new algorithm is born from a combination 
of these two algorithms. 

This newly developed algorithm has general application in both 
time-accurate and steady state flows, and also was tested extensively 
for various flow conditions, such as turbulent flows, ch emica l l y react- 
ing flows, and multiphase flows. A list of publications and doctoral 
dissertations resulting from this effort is provided in Appendix 1. 

1 CURRENT NUMERICAL METHODS 
FOR COMPRESSIBLE AND INCOM- 
PRESSIBLE FLOWS 

The objective of this study is to improve the accuracy and efficiency of 
numerical simulation of combustion ch a m ber flows for better understand- 
ing of the physics of these complicated flow conditions. To achieve this 
objective, a three year study has been initiated to develop an algorithm 
(or al g orit hms ) for incorporation of complex physical phenomena such as 
turbulence, compressibility, chemistry and multiphase effects into exist- 
ing CFD codes. This report summarizes the achievements of developing 
a unique algorithm during the last three years. Two pressure-based algo- 
rithms have been identified and extensively tested. In the following, we will 
describe the methodology and testing results. 
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The unsteady compressible Navier-Stokes (N-S) equations are a mixed 
set of hyperbolic-parabolic partial differential equations (PDE’s) while the 
unsteady incom pressible N-S equations are a mixed set of elliptic-parabolic 
PDE’s. Traditionally, different numerical techniques have been used to 
solve the N-S equations in the compressible and incompressible flow regimes. 
If the unsteady terms are dropped from the equation set for compressible 
flows, the resulting equations become a mixed set of hyperbolic-elliptic 
equations which are difficult to solve because of the differences in numeri- 
cal techniques required for hyperbolic and elliptic type equations. Conse- 
quently, nearly all successful solutions of the compressible N-S equations 
have employed the unsteady form of the governing equations. This strategy 
will also be adopted in this study for the developed algorithm to obtain a 
more general application in both time-accurate transient and steady 
state applications. With this approach, the steady state solution is ob- 
tained by marching the solution in time until convergence is achieved. 

1.1 Density-Based Methods 

To date, numerical methods utilizing the primitive variables: density, pres- 
sure and velocities (in contrast to stream function-vortidty formulation) 
for solvin g the unsteady N-S equation, can be largely classified into two 
schemes: density-velocity scheme and pressure-velocity scheme. Most den- 
sity-velocity methods have their origins in external aerodynamics problems. 
In these problems, it is natural to choose density as a primary dependent 
variable for the continuity equation, whereas the pressure is calculated from 
the equation of state. 

For most density-based methods the equations governing continuity, mo- 
mentum, energy and other scalars are solved in a coupled manner. For 
inviscid flow calculations, explicit methods (1) are often used for simplic- 
ity and storage considerations. However, the explicit methods suffer from 
the limitation of small time steps due to stability requirements, and their 
application to viscous flow problems is costly. Thus implicit methods are 
employed for most of the compressible viscous flow calculations. The most 
widely used implicit schemes for viscous flows are the methods of Beam 
and Warming [2], Briley and McDonald [3] and MacCormack [4]. Since the 
governing equations are solved in a coupled m a nne r, the characteristics of 
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the equation system are easily obtained. The fluxes at cell faces can be 
calculated by the so called flux-vector splitting techniques or by solving the 
Rieman problem. The advantage of this approach is that the jump (dis- 
continuity) conditions are satisfied, therefore good results for discontinuous 
problems, such as flows with shocks, are expected. One approach, using 
the flux-vector splitting techniques, is the Steger and Warming [5] method, 
in which the flux vector F(U) is split into two parts, F + and F~ such 
that all the eigenvalues of §£j are non-negative, and the eigenvalues 
are non-positive. The spatial derivatives of F + are backward differenced 
and that of F~ are forward differenced. The method involving solutions 
of the Rieman problem was originally proposed by Godunov (6). Since an 
exact solution of the Rieman problem is expensive and unnecessary, several 
flux- difference split schemes have been developed, for example by Colella 
[7], Dukowics [8], and Roe [9], for which the exact solution of the Rieman 
problem is replaced by an approximate solution. Since the use of central 
differencing for convective terms at high Reynolds numbers resulted in spa- 
tial oscillations in the early study, these oscillation are suppressed now by 
adding artificial damping, usually fourth-order damping terms are used by 
density-based methods. In the last few years, significant progress has been 
made in the high resolution numerical schemes based on the Total Vari- 
ation Diminishing (TVD) principle, introduced by Harten [10] to develop 
oscillation-free schemes. The TVD schemes are very robust for transient 
problems and shocks capturing [11]. Further details of these schemes can 
be found in the recent book by Hirsh [12]. 

Although the density-based coupling schemes have been successfully ap- 
plied to compressible flows, the methods have a disadvantage when in the - 
limit the incompressible flow is approached, and the linkage between pres- 
sure and density weakens in the low Mach number range [13]. In fact, Mach 
number aero constitutes a singularity in the compressible form of the equa- 
tions. Any tiny disturbances of density are enough to destroy the stability 
of solution in the low Mach number regime. Numerical experiments [14] 
confirmed the low Mach number inefficiency and instability of the density 
based methods. To avoid such a problem, a "pseudo” (artificial) compress- 
ibility term can be added into the continuity equation [15, 16]. This new 
parameter, called "pseudo”, should speed its use. As a result, the convert 
gence rate highly depends on the choice of this value. This method is not 
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efficient for unsteady simulations since many inner iterations are required 
to obtain divergence free solutions at each time step. Many existing meth- 
ods, developed specifically for incompressible flows, surmount this problem 
by treating the pressure as a primary dependent variable. These pressure- 
velocity coupling schemes are equally valid for compressible flows. 

1.2 Pressure-Based Methods 

The earliest development of primitive variable schemes based on pressure- 
velocity couplin g was the semi-implicit transient Marker-and-Cell method 
(MAC) [17] and Simplified MAC [18] by the Los Alamos group. Since there 
is no explicit governing equation for the pressure field, these methods ba- 
sically derive a working pressure equation through joint manipulation of 
the momentum and continuity equations. Existing pressure-velocity cou- 
pling methods can be divided into two categories, namely, semi-implicit 
and full implicit schemes. Because of their reliance on explicit differenc- 
ing, semi-implicit schemes have a disadvantage in time-dependent compu- 
tations, since the time-step size, necessary to maintain stability of such 
methods may drastically impair the efficiency of the algorithm especially 
when applied to the calculation of steady-state flows. Implicit methods on 
the other hand, do not suffer from time step restrictions. 

The most popular method using pressure-velocity coupling schemes for 
solving compressible flows is the SIMPLE algorithm of Patankar and 
Spalding [19] and its variants: SIMPLER by Patankar [20], SIMPLEC by 
Van Doormaal and Raithby [21], SIMPLEX by Van Doormaal and Raithby 
[22] and SIMPLEST by Sha [23]. The advantages gained by the implicit 
differencing of the SIMPLE method, which is based on a pressure correction 
procedure, are offset by the use of an iteration, which makes time depen- 
dent calculations rather expensive as iteration is required at each time step. 
The SIMPLE algorithms can be extended to handle compressible flow cal- 
culations as shown by Van Doormaal tt al. [24]. This method accounts for 
additional variations in density through an equation of state based pressure 
-density coupled correction scheme. Although applicable to a wide variety 
of flows, there are certain flow situations in which this method is inappro- 
priate and fail* to yield acceptable results. Recognized and addressed by 
Gosman and Watkins [25], these flows are the ones in which the tempera- 
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ture is strongly coupled with the pressure and velocity, such as chemically 
reacting flows. 

Another method for handling the pressure- velocity coupling of implicitly 
differenced fluid flow equations is the non-iterative PISO algorithm of Issa 
[26]. This method splits the process of the solution into a series of predictor 
and corrector steps that, at each step, a simplified set of equations in terms 
of a single unknown variable is obtained. The PISO algorithm has exhibited 
a very efficient and robust nature when applied to a variety of flows as shown 
by Benodekar [27] and Issa [28]. 

The operator splitting technique can be described briefly as follows: We 
first describe the N-S in a fully implicit way such that the N-S equation 


d , . d . . 


dp d _ 

a^ + a^ Tii + Si 


(i) 


becomes 


i {(,*., r +1 - (/>«,)"} = H (<«) - A iP " +1 + Si (2) 

where A is the discrete V» u i + \Pi +1 denote i th grid point values at (n + l)* 
time level. H( ) is a linearised convection diffusion operator. H( ) in 
principle is a time-dependent operator in which H(ti n+1 ) represents a fully 
implicit formulation while H(u n ) represents a fully explicit formulation. 
The starting point of the operator-splitting technique in this research is 
the one originally proposed by Issa [26], namely, the Pressure Implicit by 
Splitting of Operator (PISO) algorithm. In this algorithm, each marching 
time step is further divided into a sequence of predictor-corrector steps. In 
the first momentum predictor step the operator is split in such a way, 


H(u,- ) = AX + (3) 

in which A* is the diagonal part of the original matrix operator H. u* 
is the first predictor velocity value (unknown). The implicit momentum 
predictor step then can be rewritten as: 

(s-7V” u < = * ,K) ' <iP " +Si+ ^ (4) 
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Note that the pressure value used in equation (4) is the "old” value 
(time step n). Thus, the solution of (4) will not satisfy continuity at time 
step n 4- 1. 

Corrector steps are devised to derive a pressure governing equation, 
driving the velocity field to satisfy continuity. This is done by splitting the 
operator for first correcting the velocity field such that 


F«*) = A 0 u“ + Hf(u;) (5) 

thus the discretized momentum equation becomes 

(Jr7V“ r = ff ' K) “ AiP ' +Si+ ^ (6) 

Note that equation (6) now has two unknowns: «•* and p . 

By subtracting equation (4) from (6), we obtain a momentum increment 
equation 

( 7) 

By taking divergence of this equation and involving continuity con- 
straints, a pressure increment equation is obtained 

N(s - 7)“ a '} - < 8 > 

where $() is the equation of state linking pressure and density. This one 
-step corrector scheme resembles the SIMPLE algorithm, if iterations be- 
tween first predictor and first corrector were executed. To reduce this 
iteration procedure, a second corrector, based on further operator-splitting 
for second corrector velocity field, is used: 

*(«;-•) = + W) (9) 

' p m 

In this operation, the momentum and continuity then are simultane- 
ously satisfied. Summaries of the PISO procedure are listed in Table 1. 

Here, we compare the PISO algorithm to another established pressure 
based method. The newly developed FICE scheme of Hu and Wu [29] based 
on the earlier ICE scheme [30] is chosen. However, a modification of the 
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FICE scheme in parallel with the operator-splitting idea is developed as 
part of the research effort of this task in section 2. This newly developed 
algorithm is called MFICE. In the FICE scheme, the discretized momentum 
equation is written as 

(i + x) u ? +1 = ff'K) - A iP “ + S. + (10) 

The main difference between this equation and equation (4) is the split 
operator Ht(). FICE used an explicit scheme for the predictor step while 
the PISO used an implicit scheme. Another difference lies in the way the 
pressure equation is set up. Instead of deriving a pressure-correction equa- 
tion, the FICE scheme directly takes divergence of the momentum equation 
(in contrast to the momentum increment equation) and invoking the con- 
tinuity equation. This scheme is essentially a one-step predictor-corrector 
scheme and requires iteration. 


1.3 Relationship Between Pressure and Density Based' 
Methods 

Recently, the PISO algorithm has been rearranged in a vector formulation 
for direct comparison with other density based algorithms. The findings 
[31] indicate that the PISO algorithm alters the sonic speed so that the 
equations stay well conditioned in the limit of low Mach numbers. In par- 
ticular, the PISO algorithms is very closely related to the preconditioning 
algorithm developed by the Penn State group [32]. The philosophy of the 
preconditioning technique is to cause the density-based method to appear 
pressure-based at low speeds but to remain density-based at high speeds. 
Ori ginall y, preconditioning methods were used as a means of circumventing 
the disparity in the eigen values at low Mach numbers. This technique al- 
ters the tim e derivatives of the equation of motion with the acoustic speed 
scaled down to the level of the fluid velocity, such that the local CFL num- 
ber (which controls the marching time-steps) is approximately of the same 
order for viscous and invisdd terms. The extensions of the predicting tech- 
nique have recently been carried out for viscous chemical reacting flows 
involving chemical source terms. However, how preconditioning may be 
applied to improve convergence and robustness in the calculation of multi- 
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phase and turbulent flow (involving higher-order turbulent models) remains 
to be seen. 


2 NEW NUMERICAL FORMULATIONS 
FOR ALL-SPEED REGIMES BASED 
ON PRESSURE METHODS 


2.1 Formulation of PISOC and MFICE 


The basic ideas of the operator- splitting technique have been described 
in 1.2. The algorithm used was the PISO algorithm first proposed by 
Issa (26). The major difference between the PISO and the previously used 
pressure-based schemes, such as the SIMPLE-family schemes, is the use 
of momentum per unit volume as the resulting variable of the momentum 
equations, rather them velocity (which is momentum per unit mass). The 
advantages are twofold. First, the time-dependent equations give directly 
the change in a property per unit volume, whereas the SIMPLE pressure 
correction algorithm must divide the time c h a n g e by density in order to 
calculate the new velocity. Second, the change in momentum can be re- 
interpreted as a change in mass flux. This gives a linkage between pressure 
and mass flux; the mass conservation equation then only contains density in 
the time-dependent term. The pressure correction equations thus obtained 
become the momentum correction equation: 

= -(i- - p")] (li) 

thus the pressure increment equation is obtained by taking divergence of 
(11) and involving the continuity equation 




St 


](?■-?”) = 




( 12 ) 


where *(p n ,T n ) = p' Ip' 

Note that the pressure-momentum linkage equation as described in (12) 
remains essentially elliptic at all flow speeds and cannot mimic the hyper- 
bolic behavior of the system of equations, when flows are transonic and 
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supersonic. To remedy this defect, we propose to split the operator, con- 
sistent with the compressibility effect of the density correction, such that 
the first momentum corrector step becomes: 

(li - 7 ■)<>■«?* - ff 'W> - A <p' + a + 7T (is) 

the operator splitting procedure used is 

H(u;-) = H'(u;) + AA m u? (14) 

This new algorithm is named PISOC (Consistent) due to its consistent 
operator splitting procedure. Compared to the original PISO algorithm 
where 


stun = if'w) + A,^«r (15) 

the new momentum correction equation becomes 

- oX --(i-r [ Ai(p ' ” pn) + (16) 

where p! = p m — p n , and the pressure correction equation is 




(IT) 


The capability to solve compressible flow with shocks by the PISOC 
algorithm is achieved by the convection incremental pressure term, the 
second term in the left hand side of equation (17). This term properly 
takes into account the hyperbolic nature of supersonic flows. 

The PISOC algorithm is a pressure correction method (PCM) in which 
the pressure correction equations are solved. On the other hand, the 
MFICE scheme is a pressure substitute method (PSM) in which the Poisson 
pressure equations Me directly solved in place of the continuity equation. 
In the MFICE algorithm, the equation set is discretized as 
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( 18 ) 


Up"' - p") + Alw)"* 1 = o 

dt 

~ _ A 0 {pui) n + X = H/{ u? +1 ) - A<p n+l + 5" +1 (19) 

To solve the fully implicit equations (18), (19), the MFICE scheme 
employs the iterative solution procedure as shown below i 

i. k = l,(p,p,pui) fc = {p,pui) n 

ii. substituting (19) into (18) (also via equation of state p = pRT), we 
obtain the pressure equation in the form 


his* - *[(£ - = & - *{(* - a *)" 

+ #'(«?) + S? ] (20) 


which yields p k + l and thus solves for (pti*)** 1 by 

(I - A,)(^)‘ +1 = i(^r + *'(«?) + Sf (21) 

iii. if | p fc+l = p* |< e then (p,p,ptt < ) n+l = (p»P»P*«) fc=1 > g°« a to the next 
time step or else (p,p,pui) fc — (p»P»P tt *) fc+l »k = ^ + 1, goes to ii) 

This MFICE algorithm has the following characteristics: 


• The MFICE, as compared to its original FICE, holds now the operator- 
splitting (i.e, write H = A 0 + H/) scheme to enhance the convergence 
of the iteration procedure. 

• We may write H in the form H n , or H h , where the former is a linearized 
operator (depending on the n th time level results), and the latter is in 
a non-linear form, which is updated during the k — ► fc + 1 iteration. 
In MFICE Hk is used. 
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• The te chni que of introducing convection terms to the pressure equa- 
tion, to mimic the hyperbolic property when supersonic phenomena 
are involved, does not have much significance to MFICE, since the 
term dp j diva the continuity equation (1) is retained, which recovers , 
the hyperbolic property automatically. Therefore, MFICE does not 
use this technique. 

2.2 Boundary Conditions 

The implementation of boundary conditions is one of the most complex 
problems in computational fluid dynamics. The difficulties are due to var- 
ious possibilities of combining different boundary conditions in a general 
CFD problem. In addition a few special results are known about the 
mathematical representation of boundary conditions to ensure existence 
or uniqueness of the solution. For these reasons a discussion of general 
boundary conditions of CFD problems is not undertaken here. In this sec- 
tion only the most frequently encountered boundary conditions in fluid flow 
problems and the treatment, necessary to incorporate them in the discre- 
tised equations, are described. 

2.2.1 Inlet Boundary Condition 

At the inlet boundary, the values of all dependent variables are normally 
known. These are usually obtained by reference to experimental data, anal- 
ysis or estimation. These values can be substituted into the discretization 
equations for the boundary control volumes and thus nothing special needs 
to be done. Since the inlet boundary mass fluxes are generally known in in- 
compressible flows (here only the velocity component normal to the bound- 
ary is of concern), all intermediate values of Uj at the boundary, namely, 
u r* and up", are set to the given boundary value. This is equivalent to 
prescribing zero gradient boundary conditions for the pressure correction 
equation. It is readily implemented by setting the coefficient of equation 
for the boundary node to zero. The same treatment also applies for outlet, 
symmetry, and wall boundaries. The pressure at the boundary is obtained 
by linear extrapolation from inner points. 

In contrast to incompressible flows, the pressure or stagnation pres- 
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sure is often fixed at inlet boundary for compressible flows. The number 
of variables that can be specified at the inflow boundary depends on the 
number of incoming flow characteristics. For subsonic inflow, this requires 
the specification of three variables. Whereas, if the inflow is supersonic, all 
variables must be fixed. For the case of subsonic inflow there is a consider- 
able choice as to which variables should be specified. For example, both the 
velocity components and the pressure may be specified, or both the velocity 
component and temperature may be specified. In internal flows, it may be 
convenient to specify the stagnation temperature, stagnation pressure, and 
transverse component of velocity or the inlet flow angles. Other variables 
that are required at the inflow boundary are obtained by extrapolation 
from the interior. If the pressure at the boundary is specified, then p* - p" 
and p** _ p* are set to zero, which serves as the boundary condition for 
the pressure-increment equations. The normal velocities at the boundary 
are updated by the continuity equation, and the tangential velocities are 
obt aine d by linear extrapolation. If the stagnation pressure at boundary 
is specified, the velocities are obtained by the same procedure described 
above. From this velocity, the given stagnation pressure and the stagnar 
tion temperature, static pressure and temperature can be calculated. So, it 
is basically the same as the specification of a pressure boundary condition. 

2.2.2 Outlet Boundary Condition 

For i ncom pressible flows the value of the dependent variables are gener- 
ally unknown. The outlet boundary should be placed sufficiently far down- 
stream from the region of interest. As a result, any inaccuracy in estimating 
the outlet conditions will not propagate far upstream. In this study two 
conditions are used to get velocities at the outlet. The first one is that the 
velocity profile at the outlet is similar to the velocity at the first internal 
point from the outlet. The second one is that the overall continuity must 
be satisfied that is, the sum of outlet mass fluxes equals the sum of inlet 
mass fluxes. 


u° = au 1 and v° = av 1 


( 22 ) 
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( 23 ) 


rh M = ]T(u”«y„ - v°6x„) = v a(u' Syo - »'«*.) = m (n 

o 1,0 

rhjn 

a ~ T, [u I Sy 0 - v l 6x 0 ) 

1,0 

where a is & scale constant. 0 denotes values at the outlet. I denotes the 
first interior point from the outlet. When the new velocity is obtained at 
internal points, equation (24) is employed to get a, then the velocities at the 
outlet are updated by using equation (23). The other scalar variables are 
obt ain ed by setting all coefficients corresponding to the outlet boundary 
nodes to zero. It is, therefore, appropriate to evaluate outlet boundary 
values by extrapolation from upstream. 

If flow is compressible, the number of variables that can be specified 
is similar to the inlet boundary. It is equal to the number of incoming 
characteristics. For a subsonic outflow, this requires specification of one 
boundary condition. The variable that is usually fixed is the static pressure. 
It is called the back pressure. The values of other variables are obtained 
by extrapolation from the interior domain. For supersonic outflow, all the 
variables at the outlet are determined by upwind information, no boundary 
conditions should be specified, and all values of variables at the outlet are 
obtained by extrapolation. 

2.2.3 Symmetry Boundary Condition 

The symmetry boundary condition implies two contents, no flow crosses 
the symmetry line (or plane), and diffusive flux at the direction normal to 
the symmetry line is zero. The first one can be satisfied by setting the 
contravariant velocity (velocity normal to the symmetry line) to zero. The 
second one can be met by setting all the coefficients corresponding to the 
symmetry boundary points to zero. The values of variables at the symmetry 
boundary are obtained by |£ = 0, 

(*n + Vn) <t>V ~ (*<*" + = 0 ( 25 ) 
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2.2.4 Wall Boundary Condition 

For incompressible flow the no-ilip condition is applied. The velocity of the 
fluid at the wall equals the velocity of the wall. These conditions are easy to 
impose* Because the velocity at the wall is known, the implementation of 
wall boundary conditions is the same as inlet boundary conditions. There 
are two types of boundary conditions; the temperature and other scalars, 
fixed wall value and fixed wall flux. 

For compressible flow there are two possibilities of wall boundary con- 
dition, no-slip wall and slip wall. The no-slip wall is for viscous compress- 
ible flows. The implementation is the same as for the incompressible wall 
boundary. The slip wall is for inviscid compressible flows. This condition 
is imposed by setting the normal component of contravariant velocity to 
zero. The velocity at the boundary is obtained by a projection of interior 
points along the wall. 


3 EXTENSION OF THE NEW METHOD 


3.1 Turbulent Flows 

Implementation of a two-equation model (k - c), to include turbulence ef- 
fects into this combined algorithm, has been developed for subsonic flows. 
In the following, we discuss the new method and the implementation pro- 
cedure. 

The equations of fluid flows are written in general form: 


dp d . . _ 

_ + _ ( ^ )=0 

d d dp d _ 

*<'*“> + = -d^ + dr i Tii+ ‘ 


(26) 


(27) 


In the equation above p is the density, u, the velocity components, p is 
the pressure, S% are body forces and Ty are the components of a deviatoric 
stress tensor: 


Tij = P 


' du{ duj 
.dxj dxi 


39 */” 


] 


(28) 
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As mentioned in the previous section, the PISOC algorithm is a pres- 
sure corre c tion method (PCM) in which the pressure correction equations 
are solved. On the other hand, the MFICE scheme is a pressure substitute 
method (PSM) in which the Poisson pressure equation is directly solved in 
place of the continuity equation. The key issue is to include the momen- 
tum source calculations and to ensure continuity (Eq. (26)) for a prescribed 
number of corrector steps on a pressure equation. Here we seek a spht op- 
erator in a time domain such that the splitting error of the finite-difference 
form of Eq. (2) is less than the truncation error of the temporal finite 
differencing. 

The splitting procedure described above is extended for solving other 
scalar transport equations. In simulating turbulent flows, the well-known 
k-e model of turbulence requires solving transport equations for the tur- 
bulence kinetic energy, k, and its rate of dissipation e. These equations are 
strongly coupled, especially through the source terms. The splitting pro- 
cedure presented here does away with iterations, however, a non-iterative 
scheme must also be developed to deal with the other equations, such that 
the accuracy and stability of the overall scheme are preserved. It is often 
the case that the poor resolution of these scalar fields (including the species 
equation for chemical- reacting flows) undermines the integrity of the overall 
solution procedure. 

The k — e model in differential form is: 

+ = + s * (29) 

and 



where 

Si, = G — pt (31) 

and 

5« = UCtG - C,pe] (32) 

K 
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and 

G = is the generation term. The eddy diffusion coefficients 

1 J OX j ' 

p k and fit are related to eddy viscosity by 


Pt 

Pt 

fc 2 

resulting in the eddy viscosity p t = C u p— (33) 

The splitting procedure used in this study is to reconstruct the source 
terms such that 

Sk = G - (34) 

Pt 

and 

5* = C\C V — pk — c 2 cy e (35) 

pt 

In doing so, the differential equations can be split into the following: 

(f t -K.+ ^ fcn ) fc * = K, ( k ‘) + G " + e ~jT ( 36 ) 

and Predictor: 


(f t -L. + C,C^ k -) 

G P n e n 

« = Lf(e ) + Ci—Cfipk + — 

(37) 

Corrector: 



(L K. + fC-k 

\St Pt 

■)fc“ = A'/(fc*) + /x t 'G+ P ^ 

(38) 

and 



(«" L -+ c ’ c -^) 

G p n c n 

ie" = Me’) + C.-C^pk" 4- 

(39) 
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It can be shown that Jfe’* and e’ m are second order approximations in 6t 
to the original equations (25) and (26). Thus, further corrections in k and 
t are not needed. 


3.2 Chemical Reacting Flows 

The pressure-velocity procedure has also been extended to compute chem- 
ically reacting flows. In reacting flows, very large density gradients arise, 
leading to strong non-linear coupling of the equations. The incorporation 
of the species and energy equations demands significant restructuring of 
the predictor and corrector steps in the algorithm. 

The governing equations for multiple species undergoing chemical re- 
actions are the continuity, momentum, energy, and species equations. In 
generalised tensor notation they can be written: 


9 d . . A 

a, , a , , dp 

+ = ~ai 


(40) 


) 9 \ ( dui duj\ 2 Out, 

i + dxj dxi ) 3 '^dxk 


d . . d .. d 

m {ph) + d^ ipUih) = aTs 


r £] + l^ + l 


(41) 

(42) 


§i Wi) + =hi+ £; ( r f^) ’ i=1 ’ 2 ’ -’ Ar (43) 


where p is the density; u, the velocity; p the pressure; p the effective 
viscosity; h the static enthalpy; T the effective diffusion coefficient; /< the 
mass fraction of chemical species; and Ri the chemical source term. 

In addition to the equations described above, expressions are required 
for the thermodynamic quantities. In the present study, the JANNAF data 
bank [33] and the CEC76 program were incorporated for chemical equilib- 
rium and thermodynamic property estimations. In the CEC76 program, 
the minimisation of the Gibbs free energy method was used to calculate 
the composition of chemical equilibrium species. The static enthalpy and 
the fluid density are then obtained by: 
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(44) 


k=Y,fMT) 


t=l 


and 


P = P/AT|^ (45) 

where Ro is the universal gas constant; and Mwi, the molecular weight of 
species. 

For finite rate chemistry, the production rates for each of the species 
(f^) required in Eq. (43) are obtained using the laws of mass action. For 
a general homogeneous chemical reaction, which may proceed in both the 
forward and reverse directions, the stoichiometric equation can be written 
as 

f j/jjAi X>£Ai, j = l,2,...,JVn (46) 

im 1 


where Nr is the number of reaction steps. 

The law of mass continuity states that the net production of species t 

by reaction j is: 



-0 


k fi n 


V 


* S’ 

-*«ncv’ 

i=l 


(47) 


where C% is the concentration of species. The net rate of change in concen- 
tration of specie i by all reactions is found by summing the contribution 
from each reaction considered 



(48) 


In reacting flows, a coupled implicit solution procedure of a chemical 
kinetics /fluid dynamics problem would require the inversion of a complex 
system of matrices. In the present methodology, the chemical kinetics and 
the fluid dynamic solutions are decoupled in performing the integration by 
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using the operator-splitting technique [35, 36]. This procedure is embed- 
ded in the previously described predictor-corrector sequence with special 
treatment of the species equation. Using the operator representation, lets 
the governing equations be written in the following form: 

^ + C(/ i ) + D(/ i ) = ft <«) 

at 

where C() and D() are convective and diffusive operators for species /< 
respectively. 

To facilitate the splitting technique, the chemical kinetic solution only 
involves a time dependent term during the predictor step. Thus, the equa- 
tion 


^pfi = Mfuh, (50) 

at 

is integrated in a fully implicit fashion. The effective chemical source terms 
are then determined by dividing the increment of the chemical species by 
the fluid residence time. 

Subsequently, the convection and the diffusion part of the species equa- 
tions are then implicitly integrated in the correction step. The corrector 
step integrates the fluid dynamic part with the effective chemical source 
terms from the predictor step as follows: 

c(fo + wo = (^r)^ (si) 

3.3 Multiphase Flows 

The pressure velocity procedure is extended to include the dynamical equa- 
tions for spray droplets. Due to the strong coupling between two phases in 
terms of momentum, heat and mass exchange, the incorporation of the dis- 
persed droplet phase requires considerable extension of the corrector steps 
in the algorithm. These new formulations are described in this subsection. 

The general approach used for gas/droplet flows in this study is the so- 
called ’’tracking’’ approach. This approach requires formulating droplets 
dynamics in a Lagrangian frame within the continuum media for which 
an Eulerian formulation is utilized. This approach has the flexibility of 
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handling a poly-dispersed spray system and is ready for extension to dense 
spray effects, as compared to the ” two-fluid” approach. In the Eulerian- 
Lagrangian two phase approach, the governing equations for gas phase are: 

- < 52 > 


dpUi 

dt 




and for the particles: 


dp 

dxi 




dxi 



(53) 


(54) 


^ - F„ + 9i ( 55 ) 

ctt 

Since the formulation here is essentially a statistical approach, each 
computational parcel represents a large number of droplets having equal 
location, velocity, size, and temperature. The two-way coupling between 
the two phases is accounted for by the interaction terms, where 



Ui+ u'i - Vi 

* ~ r 

(56) 


S m — 51 Npm^pl dV 

psl 

(57) 

for evaporating spray 



*Wi 

ii 

N p rh tVtP (v i ) p - ^Pdr\N p fdV 

(58) 

r — - - 

in which dV denotes the computational cell, and the effective relaxation 


time r = <.//, with U = pd<%l 18 and / = C D Re P /2A. 

The goal of the present method is to build the coupling procedure of the 
two- phase interactions in a non - iterative , time - accurate frame which 
eliminates the conventional global iterations between the two-phases as de- 
picted in the Particle-Source-In-Cell (PSIC) [37] methodology. The PSIC 
method first proposed by Crowe at Washington State University in the late 
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1970’s, is still the state-of-the-science methodology and is widely used in 
industry. Our method arranges the two-phase coupling into a sequence 
of predictor-corrector steps, utilizing the operator-splitting technique, and 
thus does not require global iterations. This method is described briefly as 
follows: We seek the finite difference form of the governing equations (49) 
and (51) as follows 

E. _ 4 ,) = Hi(U p +1 ) - hP n+1 + Si + i ? +1 (59) 

and 


."+ 1 


~ v i 


Ur 1 4 - 17 / — v r i 


n+1 


St 


+ 9i 


(60) 


The superscripts n and n+1 denote time events t n and t n+l respec- 
tively. Operator A 0 and H'{ ) are constructed from the second-order up- 
wind scheme for the convection term and the central difference scheme for 
the diffusion term, respectively. The effective relaxation time T is evaluated 
at the second corrector level (**), to be defined later. The key step in our 
method is the separation of the interaction term as follows: 


FT+ 1 = -s»u? +1 + RV; ( 61 ) 

ft = -s:ur + jc w 

and so on. 

In the above splitting procedure, the terms 5” and R u are obtained by 
rearranging equation (56): 

1 Np . . 

= w E n p m p/ { st + t p) (63) 

av p= i 

and 


K = jyT, (st + r;) (vr - «i + * st). 


We now divide the predictor-corrector procedure as follows: 


(64) 
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Predictor Step 


( 


£_ 

Si 



U‘ = H' (U‘) - S iP n + Si- s:u: + RZ 

vc - v? , wr-vc 

St 9x T n 


(65) 

( 66 ) 


To strongly couple the interaction terms, we use the solved V { * and v ■ 
to evaluate r*, 5*, and i?;, such that a second approximation to the gaa 
velocity can be performed: 


(Z _ A.) Uf = H‘ {V-) - S,f + S, - S-„U- T + K (67) 

This is the essential step which eliminates the global iterations. A similar 
procedure was derived for the first corrector and the second corrector step 
involving droplet source terms in the pressure correction equations [38]. 
Turbulence effects then are added during the corrector steps. This approach 
is a unique approach , including a non-iterative feature which is consistent 
for all the physics involved. The algorithm is time accurate and requires 
no under-relaxation for all the steps involved. 


4 VALIDATION STUDIES 

We have carried out detailed linear stability analyses based on model equa- 
tions and von-Newman Fourier mode technique. Both algorithms are shown 
to be unconditionally stable with second order accuracy on time domain 
discretisation. Detailed discussions can be found in the Ph. D. disserts^ 
tions of Zhou (39) and Jiang [40]. Both algorithms have been successfully 
coded into existing CFD codes. Specifically, the PISOC algorithm is im- 
plemented in the MAST computer code for the first predictor and the re- 
maining corrector steps were established to utilize the concept of MFICE 
algorit hm . This iteration sometimes is required especially for transonic flow 
calculations. Various cases including both steady-state and time-dependent 
incompressible flows, subsonic, transonic, supersonic and hypersonic flows 
involving turbulence effects, chemistry effects as well as multiphase flows 
were studied and validated against relevant experimental data or other CFD 
results in the literature. In the following, highlights of these results will be 
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summarized. Other details involving grid set-ups, grid sensitivity studies, 
convergence histories and boundary condition set-ups can be found in the 
three Ph.D. dissertations. [39, 40, 41]. 

4.1 Driven Cavity Flow 

A square two-dimensional cavity with the top wall moving at a constant 
speed is calculated for the Reynolds number 1000 and 10,000. this case 
converges in one run in 80 CPU (CRAY/XMP) seconds with St = 1.0 
and a 51 x 51 grid system. Figure l(a&b) shows the computed streamline 
contours for Re = 1000 and 10,000 respectively. The contours are plotted 
with the same level of Ghia [42], who used a much finer grid system and 
long run time. This study demonstrates the accuracy of the current method 
and insensivity of the convergent rate due to grid refinement. 


4.2 Two-Dimensional Circular Cylinder 

The flow over a two-dimensional circular cylinder is an example of the 
external unsteady flow when the Reynolds number is not too small. A 
41 x 41 “0"-grid was algebraically generated. A fixed velocity was set at 
the outer boundary and a cyclic boundary condition was set at the front 
center I ns tantaneous stream function plot at several time steps are 
shown in Figure 2. More detailed calculations and comparisons of this case' 
can be found in [43]. 500 Steps of this calculation take 175 CPU seconds 
and 0.2 M words core memory. This study demonstrates the time-accurate 
aspect of the current algorithm. 

4.3 Backward Facing Step Flows 

Both benchmark cases of the laminar and turbulenct viscosity flows over a 
two- dimensio nal backward facing step were tested. For the l amin ar case, 
a critical case of Re = 800 is selected and the standard k - e model is 
implemented for calculating the turbulent flow case. A 81 x 51 uniform 
grid was used and the time step St = 0.5 was taken for the laminar flow 
case. 350 time steps are used to reach steady-state solutions with the 
reattachment length at 11.5 step heights which compares favorably with 
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Figure 2. Instantaneous streamline pattern for flow over a cylinder at Re = 200. 
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INS3D results. A 47 x 33 non-uniform grid is used for the turbulent flow 
case. The calculated reattachment length is 6.18 step heights; in good 
agreement with other CFD results [44]. The turbulent flow case took 60 
CPU seconds for 90 time steps (6< = 0.6). The calculated streamlines are 
shown in Figure 3 and 4 for the laminar and turbulent flow case respectively. 

This study demonstrates that the extention to include turbulence models 
does not affect the efficiency of the algorithm. 

4.4 Quasi-One-Dimensional Inviscid Flow 

The above mentioned algorithms were then implemented into the MAST 
code for one -dimensional benchmark case calculations to cover a wide 
range of Mach numbers. First, a one- dim e n si o nal symmetric nossle flow as 
calculated for subsonic, supersonic and shock flows. The calculations were 
performed on a uniform spacing of 101 grid points. Initial conditions for all 
variables were set over the whole domain as being constant and equal to the 
known analytical solution at the inlet with the exception of the exit pressure 
which was set at a fixed value. The inlet boundary condition consisted of 
extrapolating the static pressure from the interior which, together with the 
given values of total pressure and total temperature, allows an isentropic 
calculation of the other variables. At the outlet the static pressure was held 
fixed, and any unknown variable was found from the extrapolation. 

Figure 5 shows the steady state Mach distributions for subsonic, chocked 
and supersonic flows. This study demonstrates the most important features 
of the results, which are the capabilities of low to high speed flow calcu- 
lations, the sharpness of shock capturing, accuracy of shock positioning 
and speed of computation. The efficiency assessments of the algorithms 
axe s ummari sed in Table 1. Examination of Table 1 shows that the algo- 
rithms gave satisfactory shock capturing results with large CFL number 
and converged fast for all special flow regimes. 

4.5 Compressible Channel Flows Over a Circular Bump 

Validations of the PISOC algorithm for all-speed capabilities were further 
tested for two-dimensional flows. All the calculations were performed using 
a no rmalis ed time step of 0.1 for 100 time steps to reach steady state 
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STREAMLINES OF BACKWARD-STEP FLOW 
( Re=800 , uniform grid=61*51 ) 
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Figure 3. Streamline of laminar flow over backward facing step. 


STREAMLINES OF BACKWARD-FAQNG-STEP FLOW 
( RE =lE+5, K-e model non-uniform grid=l? t 33 ) 



X Direction 

Figure 4. Streamline of turbulent flow over backward facing step. 
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Figure 5. Quasi-one-dimensional in viscid, nozzel flow. 
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Figure 6. Results of bumped channel flow (Af-oe = 0.5). 
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Figure 7. Results of bumped channel flow (Af-**, = 0.675). 
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Figure 8. Results of bumped channel flow (A/-**, = 1.4). 
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Figure 9. Results of bumped channel flow (M-*, = 1.65). 
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solutions. Initial values for all variables were set constant over the whole 
domain for the flows tested, ranging from subsonic to hypersonic conditions. 

The calculated Mach number contours and surface Mach number pro- 
files are shown in Figure 6-9 for subsonic, transonic and supersonic features 
of the results showing the capabilities of low to high speed flow calculations, 
the sharpness of shock capturing, accuracy of shock positioning, and speed 
of computation. The efficiency of the algorithms is very consistent. For all 
cases tested the CPU is 0.18 ms per grid point per time step, which is a very 
fast convergence rate among the available pressure-correction methods. An 
axisymmetric invisdd flow striking a blunt body leading edge with a free 
stream mach number up to Ma. = 10 has also been tested, and the calcu- 
lated static Mach lines and the static temperature along the axisymmetric 
line are shown in Figure 10. 

4.6 Chemical Reacting Flows 

A careful validation of the present MAST code with both equilibrium chem- 
istry and finite rate chemistry has been carried out. First, an equilibrium 
chemistry example involving hypersonic viscous flow, with free stream Mach 
number of 10.0 past a two-dimensional blunt body with drcullar nose is 
tested. The Reynolds number is 8,600 and the free stream temperature is 
assumed to be 300 K. Figure 11 illustrates the iso- Mach contours of both 
the ideal gas case and the equilibrium air case. In Figure 12, the calcu- 
lated static temperatures are plotted along the flow symmetric line. The 
effects of equilibrium chemistry on the static temperature jump and shock 
location are significant in these calculations, which agree well with other 
calculations in the literature. 

The second case was the low-speed Burke-Schuman diffusion flame. The 
geometry and inlet conditions are illustrated in Figure 13. A global two- 
step Methane/Oxygen finite rate reaction model of [45] is employed for this 
case. The calculated flame temperature contours and mass fractions along 
the centerline are shown in Figure 14. The results compare very well with 
the analytical solutions reported in [47]. This study demonstrates XXX. 
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Figure 10. Results of hypersonic flow over blunt body. 
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Figure 13. Schematics of low speed diffusion flame. 
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Figure 14. Results of low speed CH A / Air diffusion flame 
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4.7 Multi-phase Flows 

For multi-phase flows, we first applied the current method for the predic- 
tion of particle dispersion controlling the turbulence with the mounting of 
different screens. The experimental set up of this case was conducted by 
Snyder and Lumley in 1971 [46]. Particle densities and sizes are chosen to 
examine the dispersion of light and small (46.5 fim diameter hollow glass), 
as well as heavy (87.0 fim solid glass) particles. Five thousand computa- 
tional particles were sampled to calculate the mean squared dispersion with 
respect to time. Comparison of the predicted and measured particle disper- 
sion are shown in Figure 15. The agreement is considered quite good. The 
current method has the flexibility of taking into account both the gravity 
(crossing trajectory effect) and the non-stokian drag law as compared to 

the continu um approach and time accuracy. 

A poly-dispersed pulsed hollow-cone spray case of practical importance 
is also chosen for the test condition listed in Table 2. The calculation starts 
at 5 mm downstream of the nozzle, and the information of particle size 
distribution and velocity distribution is directly taken from the measure- 
ments. Figure 16 shows the particle distribution plot and the gas velocity 
vectors for a 30 deg spray. With the back pressure of 1 atm, the intersection 
between the gas and the droplets trajectories is quite strong. The shape 
of the spray is no longer conical even for a very short time, and the spray 
penetration is suppressed due to the interactions of the droplets with the 
induced air flow. These flow patterns and spray shapes compared quite 
favorably with the experimental results. 

The efficiency assessment of the present numerical method is shown in 
Table 3 for the hollow cone spray case. The CPU time on a CRAY X/MP 
using the MAST code utilizing the present method for the transient spray 
calculations with 6t = 0.1ms is given. It can be seen that the amount of 
CPU time is reduced about one order of magnitude by using the present 
method. Also, the present method is rather particle number independent. 
This is due to the fact that the particles at each time step, and the source 
terms in the continuous passes, are updated for all the particles at each 
Eulerian control volume. In contrast for the TEACH/PSIC method, all 
the particles have to be tracked, and the continuous phase flow field is held 
frozen between the global interactions at each tune step. 
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Figure 15. Particle dispersion in & flow with mechanically produced turbulence. 
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radius (SMR) 
(pm) 

1 
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Table 2. Hollow cone spray 



MAST 2D 

TEACH /PSIC 


Particles 

CPU time (s) 

Particles 

CPU time (s) 

Single-orifice spray 
41x61 grid 

600 

126*9 

800 

1420 

300 time steps 

1200 

135*7 



Hollow cone spray 

31x31 grid 

400 

74*9 

800 

934 

200 time steps 

1000 

88*3 




Table 3. Efficiency assessment of numerical method. 
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The pressure-velocity algorithm developed in this study can be imple- 
mented in the framework of finite- volume, finite difference, or finite-element 
formulations. All the extensions of the current algorithm, including turbu- 
lence models, chemistry models (equilibrium and finite rate), and particle 
tracking subroutines were incorporated into the MAST code in a modular 
form. These submodel modules are stand-alone solvers and can be trans- 
ferred from one code to another with few modifications. For example, the 
chemistry module originally tested in the MAST code was transferred to 
the MFICE code which utilized staggered grids and iterative procedures 
with sucessfol applications within one month period. 

5 CONCLUSIONS AND RECOMMENDA- 
TIONS 

5.1 Summaries 

Efficient pressure-velocity coupling procedures have been developed and 
investigated for the calculation of fluid flows at all speeds. Two algo- 
rithms, PISOC and MFICE were studied, and a combined algorithm was 
implemented into an existing CFD code, MAST, with physical submod- 
els including two-equation turbulence models, equilibrium and finite-rate 
chemical reaction mechanisms and gas-droplet multiphase models. Some 
important conclusions are summarized as follows: 

1. A physics-co nsis tent pressure equation is derived to emcompass flows 
in all Mach number regimes. This pressure equation is implemented 
in the PISOC algorithm as a pressure correction method and in the 
MFICE alg orithm as a pressure substitution method. 

2. Both PISOC and MFICE algorithms are non-iteration algorithms and 
are unconditionally stable, based on linearized stability analyses. 

3. The pressure-velocity algorithm can be extended to include extra 
physical submodels, including turbulence, chemistry, and multiphase 
flows without invoking iteration procedures. 
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4. The current algorithm can be implemented on both staggered and 
non- staggered grid systems without deteriorating the computational 
efficiency. 

5. Time accuracy can be achieved within prescribed predictor corrector 
steps without invoking iterations. 

6. Ext ensions for arbitrary body-fitted coordinates do not affect the 
computational efficiency. 

5.2 Algorithm Limitations 

Although the current method is developed for general time-accurate tran- 
sient and steady state flow calculations, the time accurate aspect of the cur- 
rent method for high speed transient calculations is not firmly established. 
This aspect is strongly coupled with the spatial discretisation of the govern- 
ing equations. In addition, the finite rate chemistry procedure, currently 
implemented, is only loosely coupled with the fluid dynamics, based on the 
operator-splitting method. To accurately account for the aero-thermal phe- 
nomena, transient chemical reacting flows, such as the ignition processes, 
and the thermal property changes due to chemistry, have to be strongly cou- 
pled with the fluid dynamics through the continuity equation. With the 
present operator splitting, the current pressure-velocity algorithm cannot 
be used for time-accurate high-speed and chemical reacting flows involving 
complex chemistry. 

5.3 Recommendations 

The future work only concerns the numerical aspects of the pressure-based 
method. 

1. The pressure- velocity coupling procedure should be further explored 
to include density and temperature effects, due to complex chemistry 
and other sources, to establish strong coupling among fluid dynamics, 
thermal energy, and species to achieve time accuracy with arbitrary 
chemical kinetics. 
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2. More complex flows will result in more sparsity in time scales asso- 
ciated with various physics. Although the current method is stable 
nafng rather large time steps, adaptive time step adjustments should 
be devised, to match the solution based on optimal time steps asso- 
ciated with the relevant physics. 

3. Due to the capability of the current method to cover both incom- 
pressible and compressible flow regimes simultaneously, the pressure- 
velocity algorithm should be extended to include volume displacement 
effects due to the co-existence of liquid and gas with the same cal- 
culation grid. This will involve reformulating the momentum and 
mass fluxes across the gnd boundary and modifying the “pressure” 
equation within the grid. 

4. Currently, the solution method, such as the conjugate gradient square 
method used here, for solving the pressure equation, is still inefficient. 
Multigrid methods or multi-level error propagation methods should 
be included in the future, to speed up solving the pressure equation. 

5. For complex configurations, the pressure- velocity coupling procedure 
should be extended for implementation into multi-zone methodolo- 
gies. Coupling between calculation zones, involving different bound- 
ary conditions, should be developed to smoothly transfer information 
back-and-forth in the pressure-velocity coupling mode especially for 
unmatched grid patching. 
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Appendix 2 

Sample Inputs and Results 

A. General Input File 

The MAST family computer programs consists of a set of subroutines controlled 
by a short main program. The fundamental structure can be found in the MAST user’s 
manual version 1.0 [47]. The updated capabilities, resulting from the current study, are 
summarized in Table A.l. In the following, the updated input file descriptions and the 
handling of boundary conditions far the 2-D/Axisymmetric MAST code are described. 

The Input structure used in the MAST code utilizes blocks of optional nameti stt 
under several keywords. Each of the keywords and associated namelists are described 
here. These "keywords" are optional, i.e., skip if not needed. 

There are nine keywords built into the current MAST code. These are : GRID, 
BOUND, SOLV, PROPERTY, TURBULEN, SPRAY, REACTION, RUN and 
ENDJOB. In addition, there is one extra block called CONTROL which is used for 
identifying numeric options in the code. The CONTROL block is usually created first in 
the input file to be read into the program through Logic Unit 1. All entries are optio nal if 
certain entries require numerical values or logical values, they are entered after the entry- 
names, separated by at least one blank space. SI unit is used for required numerical values. 
Block name, entry names, the possible range of values, and a short description of the 
variable are described as follow. 

1. CONTROL Block : This block is always put at the top of the input file. 

RESTART: It activates the usage of a restart file as initial 

conditions by reading through LU=4. 

SWIRL: It activates the swirl velocity calculations. 

IMONJMON: It specifies the monitor grid point at 2-D map (IMONJMON), 

this allows the user to monitor the convergence progress for 
the selected monitor variable solved. 

The default value is (2£). 

MONU,MONV, 

MONP, MONTEMP, 

MONTK.MONTE: Only one variable can be specified as a monitor variable tracking on 

screen. Such as velocity at (x,y) location, pressure, 
temperature, turbulent kinetic energy, 
dissipation rate. The default is MONU. 
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ERRCG: 

NCGM: 

ERRM: 

INCOMP: 

COMPRES: 

OMGD: 

NCRT: 

OMGF: 


OMGT: 

PHI: 


OMGPHI: 


END: 


Convergency criterion far conjugate gradient 
Matrix solver. The default value is 0.01. 

Maximum Conjugate Gradient solver iteration number. 

The default number is 100. 

Termination criterion for steady state solutions 

using the time marching scheme. The default value is 0.0001. 

Incompressible flow calculations (=1). 

(INCOMP=0) Compressible flow options. 

1 for supersonic flows, 0 (default) for low speed 
flows. 

Number of corrector steps (dcfault=2). 

Interpolation parameter for face velocities 
1 for interpolation and 0 for averaging. 

The default value is 1. 

For supersonic t emp e ratur e field relaxation. 

Parameter in the finite difference limiter. 

1 — central, 1/2 — no name, 1/3 — 3rd order upwind, 

0 — Fromm schem, -1/2 — no name, -1 — 2nd order upwind. 
Weighting parameter for upwind scheme and other scheme. 
Numerical value of OMGPHI ranges from 0. to 1. The resulting 
scheme is weighted according to: 

PHI*( 1 -OMGPHI) + 1st order upwind scheme*OMGPHI 
End of Block input 


For example, a typical CONTROL block for running a laminar backward-facing step- 
flow using a second order upwind scheme and monitoring the convergence of U velovity 
on grid point (17,4) would have the following input block: 


CONTROL INCOMP OMEGD 0 PHI -1.0 OMGPHI 0.0 ERRCG l.E-1 ERRM l.E-4 
IMON 17 JMON 4 MONU OMGF 0.2 


A typical CONTROL input for compressible flow calculations is given in Figure A. 1. 
2. GRID Block : 


NX: 


Number of grid points in the x direction (>3). 



NY: 

Number of grid points in the y-direction (>3). 

XLEN: 

It identifies the length in SI unit of the entire physical calculation 
domain in the x direction. 

YLEN: 

The length of die entire calculation domain in y direction. 

READXY: 

It activated reading of an externally generated grid point X(IJ), 
and Y(IJ) from input file "sgridxT. The input 
format is given as: 

READ(4,9910) NX,NY 
READ(4,9920) ((X(U),I=1,NX)J=1,NY) 
READ(4,9920) ((Y(IJ), 1=1, NX) 4=1, NY) 

9910 FORMAT(2I5) 

9920 FORMAT(10E10.4) 

This input format can be modified by the user. 

UNIFORM: 

It specifies a uniform grid system. 

XDIR,YDIR: 

It specifies an X-direction grid (only 

dependent on I) or a Y-direction grid (oily dependent 

on J). Either one must be specified at the beginning of a line. 

1ST: 

The grid cell starts from 1ST index 

IEND: 

The grid cell ends with BEND index 

DST: 

The grid cell starts from DST 

DEND: 

The grid cell ends with DEND 

EXP: 

Grid space stretching factor (see the footnote *) using the power law 
formula. 

DELT: 

First grid cell size. EXP and DELT can 
be specified only once. 


* — "power law" formula. Such a grid generation can be described as follows: 

(a) Explicit "power law": 

For example, XDIR 1ST 10 IEND 25 DST 3.2 DEND 5.6 EXP 1.5, the power-law formulation 
for such a grid distribution internally calculated as: 

f t icy \ EXP 

X(I,J) = DST+ (DEND- DST) 

V 7 Iiend-istJ 

I=IST, BEND for EXP>0 

/ TCXJI^* T 

X(I, J) = DEND- (DEND- DST) 

v 7 1 ,iend-ist; 
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I=IST, IEND for EXP<0 

. If EXP = ±1, grid is uniform 

. If EXP > 1, grid is compressed close to DST and expanded close to DEND. 

. If EXP < -1, grid is expanded close to DST and compressed close to DEND. 

. EXP >1 and exp < -1 can be used to give a symmetric grid distribution. 

(b) Implicit "power-law" : 

If the user specified first grid cell size DELT instead of EXP, the grid generator 
automatically computes the stretching factor, EXP, and distributes the grid maintaining the 
first grid size DELT = X(IST +1, J) -X(IST, J) for DELT > 0 or the last grid size IDELT1 
=X(IEND, J)-X(IEND-1, J) for DELT < 0. 

For examples typical grid system used for a driven cavity flow in a square domain of 1 
meter by 1 meter with grid clustering near the wall regions and top driven lid region can be 
specifies using the power law formula with a streching factor 1.5 as: 

GRID NX 51 NY 51 

XDIR 1ST 1 IEND 26 DST 0.0 DEND 0.5 EXP 1.5 

XDIR 1ST 26 IEND 51 DST 0.5 DEND 1.0 EXP -1.5 

YDIR 1ST 1 IEND 26 DST 0.0 DEND 0.5 EXP 1.5 

YDIR 1ST 26 IEND 51 DST 0.5 DEND 1.0 EXP -1.5 


3. Property block (PROP) 


VISCOS: 

DENGAS: 

CPGAS: 

KGAS: 

PRG: 

TIN: 

PIN: 

UIN: 

VIN: 

OMEGA: 

GAMMA 

PSTAG: 

TSTAG: 


It specifies the fluid viscosity. The default value is 1. 

It specifies the fluid density. The default value is 1. 

It specifies the fluid specific heat. The default value is 1. 

It specifies the fluid thermal conductivity .The default value is 1. 

It specifies the fluid Prandtl number. The default value is 0.74. 

It specifies the flow field initial temperature. The default value is 0. 

It specifies die flow field initial pressure. The default value is 0. 

It specifies the flow field initial U velocity. The default value is 0. 

It specifies the flow field initail V velocity. The default value is 0. 

Angular momentum 

= c L 

Cy 

It specifies the flow field total pressure at the inlet The default value is 0. 

It specifies the flow field total temperature at the inlet The default value is 0. 
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Typical examples of using this input block to idendity invisid and viscous flow calculations are 
given in Figure A. 1 and A.2. 


4. Runtime block (RUN) 

DT: It specifies time step. 

DTMIN: It specifies the minimum tune step. 

DTMAX: It specifies the maximum time step. 

CFLN: It specifies the CFL number. 

The default value is 1. 

The real DT is calculated based on CFLN. If it is less than DTMIN, 
DT=DTMIN. 


NSTEP: Ma ximum time step number to be computed. 

NPR1: It specifies die screen monitor output frequency. 

NPR2: It specifies the monitor output restart file frequency. 

NEX: It specifies the example number. The user can code in 

SUBROUTINE EXAMPL. If this option is used, no other 
input data is required unless a change is desired. 

LPFAC: It specifies the monitor to output control volume face quantities. 

The default value is 0. 

LPGEO: It specifies the monitor to output the Jacobian coefficients. 

The default value is 0. 

STOP: If this keyword is specified, MAST code 

only checks the input data and no execution is done. 

For example, to run a job using time step 0.1 second for 300 time steps and to monitor the 
convergence progess every 10 time steps and save results every 100 time steps requires: 


RUN DT 0.1 NSTEP 300 NPR1 10 NPR2 100 


V Variable solution hint* (SQL Vi 

U,V: It idendify the velocity component to be solved. 

P: Pressure will be solved. 

TEMP: Temperature will be solved. 

TK, TE: Turbulence kinetic energy and dissipation rate will be solved. 

SW: Swirl velocity will be solved. (SWIRL must be activated in Selection block). 

PATC: Particle tracking is active. 


ORIGINAL s» 

OF POOR QUALITY 


58 



Examples using this input block for running laminar and turbulent flows are given in Figures 

A.1 and A.2. 

6. Turbulent block (TURBID: CT1, CT2, CMU, SME, SMK : Turbulence model constants 
CT1 = 1.44, 

CT2 = 1,92, 

CMU =0.09, 

SME =1.3, 

SMK =1.0, 

TKIN, i iilN : They specify the initial value of turbulence kinetic energy and 

the dissipation rate. The inlet TKIN is usually estimated based on 
TTQN=0.01xUlNxUIN and the measured TEIN should be used. If 
no information is known, the dissipation rate is calculated based of 
some estimated turbulence length scale SCALE. 

SCALE: If SCALE specified, TEIN is calculated 

as TEIN = CMUx TKIN 1 - 5 / SCALE 


It specifies the panicle Sauter mean radius. 

If SMR < 0, it specifies a constant particle radius. 

It specifies die built-in Chi-square droplet size distribution, 
ft specifies particle density, 
ft specifies particle te m pe ratu re. 


Particles are injected from a interval specified from grid 
(1ST, JST) to (END, JEND) 


FLOWP : ft specifies spray flow rate. 

VTNJ: ft specified particle injection velocity. 

NPTS : It specifies particle parcels per time step. 


The following example illustrate a solid-cone spray due to a point injection: 
SPRAY SMR 150.E-6 X-SQR DENPT 840 TEMP 298 
1ST 2 IEND 2 JST 2 JEND 2 NTPS 5 
VINJ 86.41 FLOWP 5.13E-3 
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x FUwmria™ block (BOUND): 


Boundary conditions must be specified in a patched form: 
1ST, IEND, JST, JEND, B.C. type, variable names. 

The variables are defined as 


1ST, IEND: The starting and finishing points in the I direction grids. 

JST, JEND: The starting and finishing points in die J direction grids. 

If IST=IEND or JST=JEND, then a line boundary condition is us ed. 
The following keywords used to specify B.C type : 

INLET: Inlet boundary condition, followed by variable names to be desc ribed: 

OUTLET: Outlet boundary condition. 

WALL: No slip wall boundary condition, followed by 

variables. Wall function of turbulence model activated. 


SYMMETRY: 

CYCLE: 

BLOCK: 


It specifies a symmetric b.c. 

It specifies a cyclic b.c. 

It specifies a blockage. Wall functions 
of turbulence model activated on block face. 


Variable name keywords (specify b.c.) 

U : x -direction velocity. 

V: y-dircction velocity. 

TK : Turbulence kinetic energy in INLET. 

IE : Turbulence dissipation rate in INLET. 

TEMP: Temperature at INLET, WALL or BLOCK. 

Q. It specifies the wall heat fluxes. The default value is 0. 

YH2, Y02, YHO, YH20, YH, YO, YH02, YH202.Y03 : 

Mass fraction at inlet. 


The default values are 0. 

Examples for this block can be illustrated by the calculations of a incompressible 
turbulent backward-facing step flow as follows: 

BOUND 

1ST 1 IEND 59 JST 1 JEND 1 WALL 
1ST 1 IEND 59 JST 35 JEND 35 WALL 

1ST 1 IEND 1 JST 16 JEND 35 INLETU 1. 

1ST 59 IEND 59 JST 1 JEND 35 OUTLET 

1ST 1 IEND 11 JST 1 JEND 15 BLOCK 
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9. 


EQLM: It specifies the equilibrium chemical reaction. 

This option requires another input file ‘chem.in’. 

The default value is 0. 

The input file “chem.in” requires die following keyword: 
ns: Number of species id be calculated 

ne: Numver of elements to be considered, 

nf: Number of species to be frozen if ’frozen’ is activated (T). 

bO: Element convervation constant, a set of *ne’ numbers. 

For example, air 1 02+3 .76N2 has 

number of element 0= 1x2=2 

number of element H=3.76x2=7.52, 

thus me conservation constants ate 2,7.52 or 1,3.76. 

chemical species 
chemical element 

frozen T to activate the frozen chemistry option 

A typical example is given in Figure A.4 

FINIR: It specifies the finite rate chemical reaction. 

The default value is 0. 

LCHEM: It specifies the the finite rate reaction model. 

The default is 0. 

LCHEM = 1 : 2-step H2 + 02 model 
LCHEM = 2 : 8-step H2 + 02 model 

INIS PE: It specifies the species be given initail values, the default value is 0. 

YH2, Y02, YHO, YH20, YH, YO, YH202, Y03, YH02: 

It specifies the initial mass fraction of species. 

The default values are 0 

RATIO: It specifies the equivalence ratio in 2-step H2+02 model 

Examples for chemistry block can be seen in Figures A.3-A.5. 

10. ENDJOB: ENDJOB identifies the end of the input file. 
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B. S ample Calculations 

In this section, five sample calculations of SSME thrust chamber flows are pi »« *n tri K 
The input files used for the calculations are shown based on the geometry and inlet 
conditions of [48] as well as the calculated iso-Mach contours and temperature contours; in 
Figures A.1-A.5. The five cases chosen are for in viscid, turbulent non-reacting rakmi— 
with the k-epsilon model, turbulent reacting flow with equilibrium chemistry, 
reacting flow with a 2-step reaction kinetics model according to [49], and turbulent 
flows with a 8-step reaction kinetics model [49], The calmlateri vacuum specific im p » 
values are also shown in the figures, which may be compared with experimental data of 
453.3 seconds. 

Finally, Table A. 1 exhibits an updated MAST code capability status from the previous 
survey date of October 1990 shown in Ref. [50]. 
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SAMPLE CASE 1 


CONTROL 


; RESTART 
GRID NX 
BOUND 


COMPRES NCRT 
OMGD 1.0 PHI 
ERRCG 1.0E-2 


3 OMGM 1 NCGM 50 


-1.0 OMGPHI 0.00 
ERRM 5.0E-6 IMON 


OMGT 0.50 OMGF 1.00 
81 JMON 10 MONU 


81 NY 41 AXISYM READXY 


1ST 1 I END 1 JST 1 JEND 41 

1ST 1 I END 81 ' JST 41 JEND 41 

1ST 1 I END 81 JST 1 JEND 1 

1ST 81 I END 81 JST 1 JEND 41 

TURBULENT TKIN 3.072 TEIN 10000 
PROPERTY VISCOS 0.0 


INLET 

SLIP Q o . 

SYMMETRY 

OUTLET 


PSTAG 20240946.90 TSTAG 3637. 
SOLV U VP TEMP TK TE 


GAMMA 1.2 GMW 10.18 


RUN DT l.E-5 DTMIN 1.0E-7 DTMAX l.E-2 

NPR1 1 NPR2 100 NEX 18 ; LPFAC 

ENDJOB 


CFLN 1.00 NSTEP 100 
; LPGEO 


CONTOUR OF MACH NUMBER 



X(M) 


CONTOUR OF TEMPERATURE 



-i ■ , - i ' i — 

-0397 0643 1643 2643 

X ( M ) 


Figure A.l Sample SSME Nozzle Flow Inputs and Results Inviscid 

ISP = 524.44 sec. 
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Y(M) 


SAMPLE CASE 2 


CONTROL COMPRES NCRT : OMGM I NCGM 



OMGD 1.0 PHI 

-1.0 OMGPHI 0. 

.00 


ERRCG 1.0E-2 

ERRM 1 . OE-5 IMON 

: RESTART 

GRID NX 

31 NY 41 

AXISYM READXY 


BOUND 

1ST 

1 I END 1 

JST 1 JEND 

41 

1ST 

1 I END 81 

JST 41 JEND 

41 

1ST 

1 I END 81 

JST 1 JEND 

1 

1ST 81 

. I END 81 

JST 1 JEND 

41 

TURBULENT 

TKIN 3.072 

TEIN 10000. 


PROPERTY 

VISCOS 9.05E-5 



OMGT 0.50 OMGF 1. 
81 JMON 10 MONU 


INLET 

WALL U 0. V 0.0 
SYMMETRY 
OUTLET 


00 


PSTAG 20240946.90 
SOLV U VP TEMP TK 
RUN DT l.E-5 DTMIN 1. OE-5 
NPR1 1 NPR2 100 

ENOJOB 


DTMAX l.E-2 CFLN 4.00 NSTEP 100 
NEX 18 ; LPFAC ; LPGEO 


CONTOUR OF MACH NUMBER 



CONTOUR CF TEMPERATURE 



X ( M ) 


Figure A. 2 Sample SSME Nozzle Flow Inputs and Results Tirrhulent- 

Non-reacting. 

ISP = 513.06 sec. 

PA nr 

OF POCR QUALITY 


SAMPLE CASE 3 


:ONTROL COMPRES NCRT 2 OMGM 1 NCGM 50 


; RESTART 
3RID NX 
BOUND 
1ST 
1ST 
1ST 


OMGD 1.0 PHI —1.0 OMGPHI 0.00 
ERRCG 1.0E-2 ERRM 1.0E-5 IMON 


AXISYM READXY 


1ST 81 I END 


I END 1 
IEND "81 
I END 81 


JST 41 


JEND 

JEND 

JEND 

JEND 


OMGT 0.50 OMGF 1.00 
81 JMON 10 MONU 


w ^ ^ w U 4 i 

TURBULENT TKIN 3.072 TEIN 10000. 
REACTION 
EQLM 

PROPERTY VISCOS 9.05E-5 

PSTAG 20240946.90 TSTAG 3637. 
SOLV U VP TEMP tk te 
RUN DT l.E-5 DTMIN 1.0E-5 DTMAX l.E-2 

-NDJOB NPR1 1 NPR2 100 NEX 18 ? LPFAC 


INLET 
WALL U 0. 
SYMMETRY 
OUTLET 


GAMMA 1. 

CFLN 4. 
; LPGEO 


V 0.0 


00 NSTEP100 



Figure A. 3 Mgjj-O NossX. Flow inputs e„ d Results ... 


ISP = 459.86 sec. 
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SAMPLE INPUT FILE 


"chem. in" 


specis # ns, elements # ne, frosen # nf 
6 2 2 
chemical element 
0 H 

chemical species 
02 H2 H20 HO O H ; 
elements conservation bo 
3 . , 8 . 
frosen ? 

F 




ORIGINAL PA&E IS 

OF POOR QUALITY 



SAMPLE CASE 4 


CONTROL 


; RESTART 
GRID NX 
BOUND 
1ST 


COMPRES NCRT 3 OMGM 1 NCGM 50 
OMGD 1.0 PHI -1.0 OMGPHI 0.00 C 
ERRCG 1.0E-2 ERRM 5.0E-6 IMON 


OMGT 0.50 OMGF 1.00 
81 JMON 10 MONU 


81 NY 41 
1 I END 


AXISYM READXY 


1 JST 

1 

JEND 41 

INLET 

.857 YHO 

0 . 

YH20 0. 


81 JST 

41 

JEND 41 

WALL U 0.0 V 0.0 

8'l JST 

1 

JEND 1 

SYMMETRY 

81 JST 

1 

JEND 41 

OUTLET 


1ST 1 I END 81 JST 41 JEND 41 

1ST 1 I END 8'l JST 1 JEND 1 

1ST 81 I END 81 JST 1 JEND 41 

TURBULENT TKIN 3.072 TEIN 10000. 
REACTION 

FINIR LCHEM 1 INISPE 
YH2 0.143 Y02 0.857 YHO 0. YH20 0. 
PROPERTY VISCOS 9 . 05E-5 

PSTAG 20240946.90 TSTAG 3637. 
SOLV U VP TEMP TK TE 

RUN DT l.E-5 DTMIN 1.0E-5 DTMAX l.E-2 

NPR1 1 NPR2 100 NEX 18 ; LPFAC 

ENDJOB 


GAMMA 1.2 GMW 10. lg 

CFLN 4.00 NSTEP 100 
: ; LPGEO 


OTNTCUR OF MACH NUMBER 





Y ( M ) Y(M) 


SAMPLE CASE 5 


CONTROL COMPRES NCRT 3 OMGM 1 NCGM 50 


; RESTART 


OMGD 1.0 PHI -1.0 OMGPHI 0.00 
ERRCG 1.0E-2 ERRM 1.0E-5 IMON 


OMGT 0.50 
81 JMON 


OMGF 1.00 
10 MONU 


3RID NX 81 
BOUND 

NY 

41 

AXISYM 

READXY 


1ST 1 

I END 

1 

JST 

1 

JEND 

41 

YH2 0.143 

Y02 

0.857 YO 

0 . 

YH 0. 

YHO 

1ST 1 

I END 

81 

JST 

41 

JEND 

41 

1ST 1 

I END 

81 

JST 

1 

JEND 

1 

1ST 81 TEND 

81 

JST 

1 

JEND 

41 


TURBULENT TKIN 3.072 TEIN 10000 
REACTION 

FINIR LCHEM 2 INISPE 
YH2 0.143 Y02 0.857 YO 0. 
PROPERTY VISCOS 9 . 05E-5 
PSTAG 20240946.90 
SOLV U VP TEMP TK TE 
RUN DT l.E-5 DTKIN 1.0E-5 DTMAX 
NPR1 1 NPR2 100 NEX 18 

ENDJOB 


INLET 
0. YH20 0. 

WALL U 0 . V 0.0 

SYMMETRY 

OUTLET 


YH 0. YHO 0.0 YH20 0.0 

TSTAG 3637. GAMMA 1.2 GMW 10.18 


l.E-2 CFLN 4.00 NSTEP 100 
; LPFAC ; LPGEO 
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CONTOUR OF TEMFERAIURE 



X ( M ) 

Figure A. 5 Sample SSME Nozzle Flow Inputs and Results Turbulent, 

8-Step Kinetics. 

ISP = 452.78 sec. 
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TYFBOFVLOW 

EQUATIONS 


PHYSICAL PROCESS 


ROCKET PROPEL! 

DISCRETIZATION 
{NUMERICAL SCSI 


MATRIX SOLVER 


1NLET/WALL i 

BOUNDARY CONX tWALL 


IPROORAM 


| CODE VALZDATDN 


PRE/POST 


PROGRAM 


ILVJu'Ji ! 11 f tfii ■ 


MISCELLANEOUS 




Sun*? D««: Oaota; 1990 Nw 

Table A.l. Updated features in the MAST code. (Dec., 1992) 
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